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Abstract. We study the evolution of dense neutral clumps lo- 
cated in the outer parts of planetary nebulae. These clumps will 
be photo-ionized by the ionizing radiation from the central star 
and change their structure in the process. The main effect of the 
ionization process is the setting up of a photo-evaporation flow 
and a shock running through the clump. Once this shock has 
moved through the entire clump it starts to accelerate because 
of the 'rocket effect'. This continues until the entire clump 
has been photo-ionized. We present an analytic model for the 
shock and accelerating phases and also the results of numeri- 
cal simulations which include detailed microphysics. We find 
a good match between the analytic description and the numer- 
ical results and use the numerical results to produce some of 
the clump's observational characteristics at different phases of 
its evolution. We compare the results with the properties of the 
fast moving low ionization knots (ansae or FLIERs) seen in a 
number of planetary nebulae. We find that the models match 
many of the kinematic and emission properties of FLIERs. 

Key words: planetary nebulae: general - Stars: AGB and post- 
AGB - ISM: jets and outflows - ISM: kinematics and dynamics 
- Hydrodynamics - Methods: numerical 



1. Introduction 

A fraction of Planetary Nebulae (PNe) show small scale struc- 
tures which differ considerably from their surroundings in ei- 
ther line ratios or radial velocity or both. When these appear 
in pairs in along the symmetry axis of the nebulae they are 
often referred to as 'ansae', or when there is a clear velocity 
difference, by the acronym of FLIERs (Fast Low Ionization 
Emission Regions), first introduced by Balick et al. (1993). 
The defining observational properties of FLIERs are that they 
should be bright in low ionization lines, located along the ma- 
jor axis of the PNe, and have a higher velocity than their envi- 
ronment. There are many cases which fulfill the first two cri- 
teria but for which no kinematic data is available, leaving their 



st atus as FLIERs uncer tain. Examples of these can be found 
in |Corradi et al. (1996)j who took ratios of Ha+[N II] to [O 



III] to bring out many low ionization regions in PNe. Often 
the ansae are only marginally resolved in ground based im- 
ages; in some cases they are clearly elongated into 'tails'. In 
some other cases o ne observes a wh o le series of ansae, as for 
example in Fg 1 ( Lopez et al. 1993 , palmer et al. 1996 ). See 
Mellema ( 1996)1 and |L6pez (1997)| for reviews on what appear 
to be collimated outflows from PNe. 

Among the best studied FLIERs are the ones in NGC 3242, 
NGC 6543, NGC 6826, NGC 7009, and NGC 7662, and these 
are the ones that have also been selected for imaging with 
the HST ( |Harrington 1995| ; |Balick et al. 1997[ ). Ground based 
studies of the kinematics of the FLIERs in these PNe can be 
found in [Balick et al. (1987)| and |Miranda & Solf (1992)| . The 
radial velocities of these FLIERs are typically around 30 to 
50 km s _1 . IBalick et al. (1993)1 and iBalick et all (1994) stud- 



ied their spectroscopic properties in some detail. They found 
that the derived densities and temperatures are not very dif- 
ferent from those of the environment, a surprising result given 
their conspicuous appearance in the low ionization lines. The 
higher ionization species are found at the side facing the star, 
indicating that photo-ionization determines the ionization strat- 
ification. Because of the observed velocities one would expect 
a bow shock (and higher ionization species associated with it) 
at the end facing away from the star, but this is not observed. 

The HST results^ (a preview of which can be found in 
Weinberger & Kerber 1997) show a considerable amount of 
detail and the analysis has not yet been completed. The images 
seem to show little arcs either pointing towards or away from 
the star, and more diffuse regions also pointing both towards 
and away from the star. 

The actual configuration of the FLIERs is not at all clear. 
The groundbased observations suggested fast moving clumps 
being photo-ionized from behind, but the morphologies seen in 
the HST images do not seem to support this. The tails pointing 
away from the stars (very clearly seen in NGC 3242), suggest 
rather the reverse: a clump being hit by a wind from behind. 
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However there is no other evidence for such a wind. At this 
point it seems hard to propose a definite model for the FLIERs. 

One of the surprising properties of the FLIERs is their 
brightness in [N II]. The [N II]A6583 line is normally as strong 
as or stronger than the nearby Ha line. This is usually inter- 
preted as being due to a high nitrogen abundance, which im- 
plies that the FLIERs were formed from stellar material en- 
riched through the CNO-cycle. 

To study the structure of FLIERs we developed a hydrody- 
namics code which includes the most relevant micro-physical 
processes, such as radiative cooling and photo-ionization. But 
as was stated above it is not clear what to choose for initial con- 
ditions. There appear to be two cases which could be used as a 
FLIER model. Both involve a dense clump moving away from 
the central star, but in the one case the clump is moving faster 
than its environment, and in the other slower. The first case is 
the 'classical' FLIER model, but the new HST results support 
the second model. 

What both cases have in common is the gradual photo- 
ionization of the clump. We therefore decided to study this ef- 
fect first, for the moment neglecting the effect of any velocity 
difference with the environment. In fact, the photo-ionization 
will cause an acceleration of the clump due to the so called 
rocket effect (Oort & Spitzei 1995). If the rocket effect is very 



efficient, it might not even be necessary to invoke a wind; the 
clump would reach FLIER-like velocities within its life time 
(~ 1000 years, see Balick et al. 1987). 

We therefore study here the evolution of a dense clump be- 
ing photo-ionized by the stellar radiation. In Sect. 2 we present 
an analytic model for the shape and evolution of such a clump. 
In Sect. 3 we describe the numerical method that we used to 
compute the models shown in Sect. 4. In Sect. 5 we compare 
the numerical results to the analytic ones, and in Sect. 6 we 
compare them to the observations. In Sect. 7 we summarize the 
conclusions of the paper. 

2. Analytic model 

2.1. General description 

The evolution of a dense clump being photo-ionized con- 
sists of two phases, the collapse or implosion phase and 
the cometary phase. These processes have been studied by 
many authors, starting with Oort & Spitzer (1955). The three 
most recent papers are those of Bertoldi (1989) on the col- 
lapse phase, Bertoldi & McKee (1990) on the cometary phase 
and Lefloch & Lazareff (1994) on both. The introduction in 



Bertoldi (1989) contains a good review of all the previous 



work. 

In the collapse phase an ionization front of type D forms, 
preceded by a shock wave which travels through the clump, 
gradually compressing, accelerating and heating it. At the ion- 
ization front, material is ionized and thrown out in a photo- 
evaporation flow. Once the shock has moved through the entire 
clump, it enters the cometary phase in which the whole clump 
starts to accelerate due to the rocket effect caused by the photo- 
evaporation flow. 



2.2. The collapse phase 

Let us consider the initial acceleration of a spherical clump by 
the interaction with the impinging ionizing radiation field. The 
clump has an initial mass Mi, density p u radius P; and sound 
speed q. Let M s be the mass of the shocked material and v its 
velocity (see Fig. 1). The equation of motion of this mass is 



dv 

~dt 



PqAi - p[V s A s v . 



and the equation of mass conservation is 



dM s 
dt 



-mFAi + PiV s A s , 



(1) 



(2) 



where v s is the shock velocity, and A\ and ^4s are the effective 
areas of the ionization front and the shock front (respectively), 
shown schematically in Fig. |] Po 1S the pressure driving the 
shock, which results from a D-critical ionization front. Po = 
2Fmci (see Sect. 2.3), in which F is the ionizing photon flux 
reaching the ionization front and C\ is the isothermal sound 
speed of the ionized material. We assume P to be constant. 
The average mass per atom or ion is taken to be m « 1.3 toh 
(where mn is the mass of hydrogen). If Ai w A s , Eq.[| has the 
trivial solution v = vq = constant, with quantities given by 



P = PiV s V , 

PiV s = po (v B - v ) , 

Pi (q 2 + v s 2 ) = po [c 2 + {v s ~ v ) 2 ] 



(3) 
(4) 
(5) 



where po is the density of the shocked material, and Co its sound 
speed. This sound speed is fixed by the radiative losses in the 
post-shock region. If we now assume that q <C v s , Eqs. ^ and 
| give 



v s sa v H , 



P = Pi 1 + — 

co z 



(6) 



(7) 



Substituting Eq. ^ into || one finally obtains the (constant) ve- 
locity vo of the shocked material, which is given by 



, -Po 2 

Vo = C 

Pi 



1/2 



(8) 



so that the position of the clump as a function of time is Xjp — 
Vq t. The shock goes through the diameter of the clump in a 
time 



to = 



2 Pi 



(9) 



after which the cometary phase starts (described in the next 
subsection). We now show that by time t a substantial amount 
of the clump material has already been photoionized. 
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Fig. 1. Schematic diagram showing the interaction of a dense 
clump of initial radius i?i with an impinging ionizing pho- 
ton flux Fq. The back-reaction of the photo-evaporation flow 
(driven out at a sonic velocity ci) drives a shock (of velocity v s ) 
into the clump. The shocked clump material moves away from 
the photon source at a velocity v. The position of the ionization 
front along the symmetry axis is xip. The ionization front is 
indicated with a thick line, the shock front with a dotted line 

Using the A\ k, As approximation in Eq. ^, one can con- 
clude that the ratio between the neutral clump mass at the end 
of this phase, Mq, and its initial mass is 

M _ mF 
Mi piv s ' 

The mass M of the clump can be obtained from the differ- 
ential equation 



(10) 



dM 

= -mF Ai . 

dt 



(ID 



In order to integrate this equation, let us consider the length I 
associated with the area A\ : I oc Ai 1 ^ 2 . Then, for a uniform 
density clump we have M oc I 3 , and therefore we can write 



(12) 



where we assume that f3 is constant with time. However, this 
assumption of a constant (3 is only strictly true for a clump 
of time-independent average density, which is not the case as 
the clump has both undisturbed and shocked components (see 
Fig. [j]) of different densities, and the proportion of mass in 
these two components evolves with time. However, as we will 
see when comparing the analytic model with the numerical re- 
sults in Sect. 5, this approximation appears to give a reasonable 
description of the evolution of the clump in the collapse phase. 
Using Eq. we can integrate |TT|to obtain 

3 



M{t) = Mi 1-/3 



m F 
'3Mi 1/3 



where we have taken M = Mi for t = as the initial condition. 

To find the constant j3, we use the condition M = Mq 
(given by Eq. HQ) for t = to- The result is 



3Mi 1/3 
mFto 



1 



Mi 



1/3' 



(14) 



Substituting Eq. |14j into |13|, we finally obtain the mass of the 
clump as a function of time 

^ 3 



M(t) =M; 1- 



1 - 



Mo 
M 



1/3' 



(15) 



where the term Mo/M; can of course be substituted by the right 
hand side of Eq. [l0| The flux F during the collapse phase is 
approximately given by Eq. [34] (see below) with = Ri- 

2.3. The cometary phase 

In this section, we present a description of the 'cometary 
phase', which begins when the shock has moved through the 
entire clump. At this time, to (Eq. ||), the shocked clump has a 
mass Mo (Eq. [l(J), and is moving at a velocity vq (Eq. fy. 

Let us assume that the shocked clump has a uniform tem- 
perature T (and therefore a uniform isothermal sound speed 
c ) and that it is subject to a uniform acceleration a. The pres- 
sure and density stratifications then follow from the condition 
of hydrostatic equilibrium 

— — = -ap , (16) 
dx 

which gives pressure and density stratifications of the form 

P Q {x) =p (x)c 2 = Po(0)e-*/\ (17) 

where Po(0) is the pressure at the tip of the clump (at x = 0), 
and h is the scale height, which is given by 

(18) 



h= C ^ 



The x-axis points away from the head of the clump in a di- 
rection parallel to the ionizing photon flux (see the schematic 
diagram of Fig. |lj). 

Let us now assume that the clump is embedded in an im- 
pinging ionizing radiation field, so that it is enveloped by an 
ionization front which moves into the clump. Let uq and u\ 
be the velocities of the neutral and ionized flows with respect 
to the ionization front. The mass and momentum conservation 
equations then give us 



PlU\ = p u , 

Pi + piu 2 = P + po u 2 



(19) 
(20) 



where Pi and pi are the pressure and density in the ionized 
material, and Po and po are the pressure and density on the 
neutral side of the front, which are given by Eq. |l7|. From Eqs. 
and On, we then obtain 



(13) P (x)=P 1 



P\Ul 



1 



III 



(21) 



We define ci to be the isothermal sound speed in the ionized 
material, and assume that ui — c%, as would be the case for a 
D-critical ionization front. Then Pi = p\Ui 2 , and therefore 



Ul 



1 



m 2 - 


u 1 






^ 2 





< 1. 



(22) 
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Eq. [2l] can then be written as 

P (x) « 2pu ci , 



(23) 



where we have also used Eq. [T?J 

Let us now consider a point on the surface of the clump, 
as shown in the schematic diagram of Fig. |l[ Considering that 
each ionizing photon impinging on the ionization front results 
in a photoionization of a further atom of the clump, it follows 
that 



Fm sin a 



(24) 



where F is the (uniform) impinging photon flux, m is the aver- 
age mass per atom of the clump (« 1.3 itih) and a is the angle 
defined in Fig. [j]. In this we ignore the differences in the ion- 
ization structure of hydrogen and helium. Substituting Eq. ^4] 
into E3l we obtain 



Pq (x) — 2 F c\ m sin a , 
and from Eq. |l7| 

P {0) = 2F Cl m, 



(25) 



(26) 



leading to the following differential equation for r, the cylin- 
drical radius of the clump at a position x, 



dr 
dx 



(e 2x ' h - l) 



-1/2 



Integrating Eq. |27| with the boundary condition r = for x = 
0, one obtains 



h arctan 



^ e 2x/h _ I 



(28) 



which gives the shape of the clump. This solution for the 
shape of a photo-evaporating clump is very similar to the 
solution for a clump accelerating due to a passing wind 
( de Young & Axford 1967| ), normally referred to as a plasmon. 
We will therefore also refer to the photo-evaporating clump in 
the cometary phase as a plasmon. 

The net force T x in the x-direction and the mass M of the 
clump follow from the integrations 



TTh/2 



2nr Po(x) dr , 



M = I irr pq (x) dx . 
/o 



Using Eqs. |l7|and |28[ one can integrate Eqs. 
tain 



r, = li--) TT 2 p (o)h 2 , 



M = I 1 

7T 



>^o(0) 



f'o- 



(29) 
(30) 

| and |(] to ob- 

(31) 
(32) 



which enable us to recover Eq. [L8] by taking the ratio a = 
T x jM = Co 2 /h. Equation [32] gives the mass of the clump in 
terms of the scale height (or viceversa). 



The ionizing photon flux F reaching the ionization front 
is in general smaller than the flux Fq delivered by the star at 
the position of the clump. The flow produced at the ionization 
front, streaming out of the clump towards the star absorbs some 
of the impinging photons. The difference between Fq and F 
can be expressed as ( |Spitzer 1978 ), 

a R F 2 R 1F 



Fn-F 



3ci 2 



(33) 



where ctR is the recombination coefficient, and i?ip is the ra- 
dius of the curved ionization front. Although Eq. [33] has an ex- 
act solution for F in terms of Fq and in order to proceed 
analytically with the equations of motion for the clump it is 
necessary to use the approximate solution 



F: 



F Q 



1 



orFoRif, 



1/2 



(34) 



which deviates from the exact solution of Eq. B3[by less than 
15%. For the radius of the ionizing front in Eq. [34], we will use 
the radius of curvature of the tip of the clump, which is equal 
to the scale height h. 

The initial mass Mo and the corresponding scale height ho 
are related through Eq. [32], which gives 

h 3 

^ IV / ' r I / " I III I =^^^^= 



(27) Mo=(l 



2\ 

- —(2 Foam) 

7T J Cq z 



1 



3ci 2 



1/2 



(35) 



where we have also used Eqs. ^ and gj. 

The rate of change of the mass of the plasmon is related to 
the impinging ionizing photon flux through 

dM ft \2 



= — b m it —h 

dt \2 J 



(36) 



where irh/2 is the outer radius of the projection of the plasmon 
onto a plane perpendicular to the symmetry axis (cf. Eq. prj ). 
Equating Eq. [36] to the time derivative of Eq. ^2] and using 
Eq. |34], one finds a differential equation for the time evolution 
of the scale height h 



dh 
dt 



6{l + yh) 
6 + hr\h 



in terms of the dimensionless variables h 
where 



V 



v*ho 

CO 2 

aRFphp 
3ci 2 



- ci , 

7T .' 



(37) 

h/ho and t = t/t*, 

(38) 
(39) 
(40) 



Equation |3_7j can be integrated with initial conditions h 

i = t n (= t 1 't*) to obtain 



1 - h 



1 

6r) 



In 



1 



n 



l + r/h 



t-tn 



lfor 



(41) 
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It follows from Eq. ^ijthat the clump becomes fully ionized 
(that is, h = 0) at a time 



tm=t + l + ^- 111(1+7?) 

6 bij 



(42) 



We find that Eq. can be inverted for h in an approximate 
way, obtaining 



in 



t 



tr 



to 



The mass of the clump as a function of time then follows from 
Eq. p2|, which gives 



\1 + Vhj 
with h given 

byEq.||. 
The equation of motion of the clump is 

dv 



(44) 



(45) 



where v is the velocity of the clump along the a;-axis. Using 
Eqs.[l|and|43|,Eq.|4| can be written as 



dv 
di 



1 

h 



(46) 



where v = v/v*. With h given byEq. ||Eq.|| can be inte- 
grated to obtain 



v = v Q + (i m - t ) In 



(47) 



Note that this expression means that the clump reaches infinite 
velocity at t = t m , but that this limit is approached logarithmi- 
cally. The dimensionless position of the clump x = x/(v*t*) 
follows from the time integration of Eq. which gives 



x = v t + (i m - t ) 



(tm - t) In 



{t - t ) + 



to 



Taking the limit of Eq. |48| for t 
where 

x m = v t m + [t m - t ) 



(48) 

i m we find that x — > x m , 
(49) 



is the position of the clump at the time that it becomes fully 
ionized. The velocity averaged over the lifetime of the neutral 
clump is 



< v >= v + 



to) 



(50) 



Equations ^] and |l5| give the evolution of the clump in 
the collapse phase (in which the shock has not yet travelled 
completely through the clump). The fully shocked, cometary 
regime is described by Eqs. fPu O and KM We find that for the 



cometary regime it is also possible to derive analytically the 
shape of the plasmon (Eq. ^8|). In Sect. 5 we carry out a de- 
tailed comparison between this analytic description and a full 
gasdynamic simulation of the flow. 

Our description of the cometary regime is quite similar to 



the analytic solution derived by Bertoldi & McKee (1990) The 



model described by these authors differs from the present one 
in the following ways 



(43) _ 



while the present model considers that the incident radia- 
tion field has a planar geometry, the dilution due to a spher- 
ical divergence was included by Bertoldi & McKee (1990) 



- while Bertoldi & McKee (1990) consider the effects of the 



photo-evaporation flow on the impinging radiative field 
only in the high absorption limit (where F oc i?ip -1 ^ 2 ), our 
model includes the transition to the low absorption F w Fo 
regime that occurs for low values of Rip (see Eqs. [33] and 

- the cases in which magnetic pressure and the thickness of 
the ionization front are important were also considered by 



Bertoldi & McKee (1990) 



The first two points are appropriate for the conditions in PNe, 
and necessary to allow quantitative comparisons with the re- 
sults from the numerical simulation, which are described in the 
following sections. 

3. Numerical method 

To numerically study the evolution of a photo-evaporating 



clump we used the code described in Raga et al. (1997a) 



adding photo-ionization to the processes included in the rou- 
tines dealing with the radiation and ionization physics. This 
code (called CORAL) solves the Euler equations by using a 
Van Leer flux -splitting algorith m combined with an adaptive 
grid method ( Raga et al. 1995 ). The radiation physics part is 
dealt with using operator splitting, meaning that it is solved for 
separately from the hydrodynamics each time step. 

We follow the non-equilibrium evolution of the ionic abun- 
dances of H I— II, He I— III, C II-VI, N I-VI, O I-VI, Ne I-VI, 
S II-VI. We assume that C I and S I do not occur due to the 
abundance of photons below 13.6 eV. Some of the ionic abun- 
dances are derived by normalizing to the total abundance of the 
gi ven element (He HI , C II, N I, O I, Ne I, S II). Compared 
to Raga et al. (1997a) we have added nitrogen to the list of ele- 



ments included since the nitrogen lines are so prominent in the 
FLIERs. 

The ionic abundances are determined by the recombi- 
nation rates due to radiative, dielectronic and charge ex- 
change recombinations, and the ionization rates due to col- 
lisional ionization, charge exchange, and photo-ionization. 
For the metals all of these except the photo-ionization rate 
are found from two-dimensional tables in temperature and 
electron density, calculated using the atomic data speci- 
fied in feagaetal. (1997a) For nitrogen we used the at omic 
data from Lundqvist & Fransson (1996) , Nahar (1995) , and 
|Nahar (1996) The treatment of the photo-ionization is detailed 



6 



Garrelt Mellema et al.: Photo-evaporation of clumps in Planetary Nebulae 



below. For H and He we use analytic fits to the radiative recom- 
bination and collisional ionization rates, and do not consider 
charge exchange. 

From the rates we first calculate the new ionic abundances 
of hydrogen and helium, and from these the electron density. 
This is then used to calculate the new ionic abundances for the 
metals iterating a series of two-ionization state solutions. 

Using the updated ionic abundances we calculate the heat- 
ing rate due to photo-ionization of H I, and He I— II, and the 
cooling rate due to line cooling of H I — II, C II-IV, N I-V, 
O I-V, Ne II- V. Similarly to the ionization and recombina- 
tion rates, the latter are stored in two-dimensional tables in 
temperature and density. The tables for H, C, O, and Ne can 



be found in 
piled from 


Raga et al. (1997a). Tables for nitrogen were com- 


Lundqvist & Fransson (1996), 


Brage et al. (1995), 


Peng &Pradhan (1995) 


and 


Fleming et al. (1996) 


All tables 



are available from the authors upon request. The heating and 
cooling rates are used to calculate a new temperature. This 
whole procedure is iterated several times in order to find the 
average value of the temperature during the time step. 

3.1. Photo-ionization processes 

The photo-ionization rates are found from the integrals over 
the stellar spectrum weighted with the photo-ionization cross 
section of the species under consideration. The heating rates 
are similar but have the extra factor of the excess photon 
energy over the ionization threshold value v t . The treatment 
of these integrals is similar to the approach described in 



Frank & Mellema (1994) We basically calculate the values of 
the integrals for a range of optical depths for three frequency in- 
tervals: z/ t (H I)-f t (He I), i/ t (He I)-f t (He II), and v t (Ue II)- 
oo. In principle one has to deal with two optical depths (r(H I) 
and r(He I)) in interval 2, and three (those and r(He II)) in 
interval 3. But it is possible to represent these sums with one 
optical depth by using the fo llowing approximation taken from 
Tenorio-Tagle et al. (19831 



t 2 (» 



73(f) 



[T t (HI)(0.63i/ t (HI)) 1 - 7 
TtCHelJi/^CHel)]!/- 



1.7 



(51) 



[r t (HlK 2 - 8 (HI) + 
Tt(HeI)(1.51^t(HeI)) 2 - 8 
7i (He UK 2 ' 8 (He II)] 



2.8 



(52) 



in which the r t and the v t are the values at the ionization thresh- 
old of each interval. These formulae assume that the cross- 
sections behave as powerlaws with index 2.8 for H I and He II. 
For He I we use an index of 1.7 which is accurate between 
24.6 eV and 45 eV, and at 5 4.4 eV gives only a ~ 15% overes- 
timate of the cross section (Samson et al. 1994). 



Storing the values for the three frequency ranges for 40 
optical depths ranging logarithmically from 10~ 4 to 10 4 , 
we only need to calculate the integrals once and afterwards 
interpolate in the tables. By using the approximation of 



Tenorio-Tagle et al. (1983) we correctly describe the harden- 
ing of the spectrum close to an ionization front and the effects 
of the presence of He, and thus improve on the approach de- 



scribed in Raga et al. (1997b), without having to recalculate the 
integrals every time step. The necessary photo-ionization cross- 
sections were taken from psterbrock (1989) and Cox (1970) 



4. Numerical results 

We ran a simulation of the evolution of a clump using the 
method described above. As parameters we chose an initial 
clump density of 5000 cm -3 , temperature 100 K, distance to 
the star 10 18 cm= 0.32 pc and diameter 2 x 10 16 cm. With these 
parameters the initial mass of the clump is 2.3 x 10~ 5 Mq. 
The clump is initially neutral, except for the C II and S II (see 
Sect. 3). The environment has a density of 100 cm~ 3 , tem- 
perature 10 4 K, and is initially fully ionized. The star has a 
black body spectrum with an effective temperature of 50 000 K, 
and luminosity 7000 Lq, typical values for the central stars 
of young PNe. For abundances we used the PN values in Ta- 
ble 5.13 from Osterbrock (1989) The size of the adapative grid 
at maximum resolution is 1024 x 256. 

We followed the evolution of the clump for 950 years af- 
ter which it has been completely evaporated. We find that the 
initial collapse phase lasts about 300 years, and that the rest of 
the time the clump can be described as an accelerating plas- 
mon. Here we show one snapshot from the collapse phase 
(t = 250 years, Fig. 2) and two from the cometary phase 
(t = 500 and 750 years, Figs. 3 and 4). These figures show 
the density and velocity field for these times. The lower plot 
also shows the extent of the clump material (mainly due to the 
photo-evaporation) and the areas which are neutral. The ve- 
locity arrows are only plotted if they correspond to velocities 
higher than 5 km s _1 . 

At t = 250 years, the clump which originally extended 
from i = 2x 10 16 to x = 4 x 10 16 cm, has evaporated about 
half its material. The ionization front is at 3 x 10 16 cm and the 
shock front at 3.6 x 10 16 cm. One sees the photo-evaporation 
flow at the back of the clump; th flow's density drops as it ex- 
pands and accelerates away from the clump, and there is a dou- 
ble shock structure separating it from the environment material. 
The area in the shadow of the clump is recombining and cool- 
ing, which means that it is no longer in pressure equilibrium 
with the rest of the environment material. It is slowly being 
compressed and heated; this is more apparent at later times. 

The compression through the Shockwave is about a factor 

5, leading to densities of about 25 000 cm -3 in the shocked 
part of the clump. The temperature in that part is a little under 
10 4 K, not enough to produce appreciable collisional ioniza- 
tion, and the material is therefore almost completely neutral. It 
is accelerated to about 10 km s _1 . 

At t = 500 years, the shock wave has travelled through 
the entire clump. The shadow region has become more com- 
pressed, and there is a small part at the front of the plasmon 
which actually has formed from environmental shadow ma- 
terial. The velocity arrows show that the clump has started 
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Fig. 2. Density and velocity at t = 250 years. The velocity arrows are only plotted for velocities larger than 5 km s _1 ; in the 
velocity plot the solid line deliminates the extent of the original clump material and the dashed line is the 50% H ionization 
contour. The dotted demicircle in the top figure indicates the original position and size of the clump. The contours in the density 
plot are factors y2 apart; the lowest contour is 0. 1 cm~ 3 



to accelerate, as is also witnessed by the fact that the photo- 
evaporation flow on the top of the plasmon leaves it at an an- 
gle. In the rest frame of the clump the evaporation flow is per- 
pendicular to the surface. More details about this phase can be 
found in the next section where we compare it to the analytic 
solution. The plasmon density ranges from about 2 x 10 4 cm -3 
near the ionization front to 1500 cm -3 at the far end; temper- 
atures vary from 10 4 K to 2000 K. Velocities vary from 13 to 
20 km s _1 . Due to the interaction with the environment ma- 
terial the shock in the photo-evaporation flow has acquired a 
conical shape. 

At t — 750 years, the clump is close to complete evapo- 
ration. The neutral region of the flow has become much more 
extended due to the recombination of the environment material 
in the shadow of the clump. Otherwise the situation resembles 
the one at t = 500 years. At the tip the density is still high 
(1.5 x 10 4 cm -3 ), but in most of the plasmon it is a factor 10 
lower than that. The temperature at the tip is about 6000 K, in 
the rest of the plasmon it varies between 2000 and 3000 K. The 
velocity throughout the plasmon is about 20 km s -1 . 

After this the clump quickly becomes thinner and thinner, 
until it completely disappears (between 900 and 950 years). 
The neutral former shadow region stays on for about 50 years 
more, after which it becomes photo-ionized again, and all the 
material on the grid is ionized. This accounts for the last data 
point in Fig. 5. 



We find that the photo-evaporation flow has a temperature 
of ~ 3000-5000 K at the base. Within the ionization front, the 
temperature has a sharp peak of more than 10 4 K, a result of 
the hardening of the spectrum within the front. 

5. Comparison between the analytic and numerical models 

To compare the numerical results from the previous section 
with the analytic description from Sect. 2, we determined the 
position and mass of the clump as a function of time. The po- 
sition of the clump is defined as the position on the symmetry 
axis where n(HII)/n(H) = 0.5, and the mass is all of the 
clump material with n(H II)/n(H) < 0.5. The squares plotted 
in Fig |] are the values we derived. After 900 years there is no 
neutral clump material left. 

In order to carry out a comparison with the analytic model, 

~and|48 
and 44 



in Fig. 5 we also plot the results obtained from Eqs 
(giving the position as a function of time) and Eqs. 115 
(giving the mass of the clump as a function of time). 

The free parameters of the analytic model are the initial 
density and radius of the clump (which we set equal to the 
corresponding values of the numerical simulation), and the 
sound speeds Co (of the shocked clump) and c\ (of the photo- 
evaporation flow). As discussed above, from the numerical sim- 
ulation we see that the shocked clump has a temperature of 
~ 6000 K, and that the base of the photo-evaporation flow 
has a temperature of ~ 3000 — 5000 K. From these temper- 
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Fig. 5. Position and mass of the clump against time. Squares 
are values from the numerical simulation, the solid lines are 
the results of the analytic theory 



atures, we derive values Co ~ 6.2 km s _1 and c% w 8.4 km s _1 
for the shocked clump and the ionized flow, respectively. For 
these parameters, we obtain to = 415 years, t m = 985 years, 
ho/Ry = 0.344, and an average velocity < v >— 26 km s _1 
over the lifetime of the neutral clump. 

We also compared the shape of the plasmon with the ana- 
lytic expression given by Eq. 28 The left-hand side of Fig. 6 



shows contour plots for the pressure, temperature and H ioniza- 
tion fraction at t — 500 years with a fit according to Eq. fus- 
ing h = 2.7 x 10 15 cm. The steep gradient seen at 5.5 x 10 



cm 



corresponds to the transition from the clump to the shadow be- 
hind it. The fit for the shape is very good, especially when one 
considers the complex structure seen inside the plasmon. To 
check whether the derived scale height corresponds indeed to 
the real one, we have plotted a cut of pressure along the sym- 
metry axis and overplotted the gradient expected from a value 
of h = 2.7 x 10 15 cm. This is shown in right-hand side of 
Fig. 6, together with cuts for the temperature and H ionization 
fraction. One sees that logarithmic pressure gradient is not con- 
stant but on average corresponds quite well to the value used to 
the make the fit to the shape of the plasmon. 

Another method to derive the scale height is to use Eq. [|3| 
Using t = 500 years and the values of t , t m and h de- 
rived above, we obtain a value for the scale height h = 
2.9 x 10 15 cm, in excellent agreement with the value obtained 
from the fit to the shape of the plasmon. 

Given the very simple nature of the analytic model and 
the many details included in the numerical model, the match 
between the two is surprisingly good. Both the temporal be- 
haviour (Fig. 5) and the shape of the plasmon (Fig. 6) found in 
the numerical simulation closely match the analytic ones. This 
can be seen both as a test of the numerical method and a confir- 
mation that the assumptions going into the analytic models are 
largely valid. The match is even more surprising given the poor 
match found between the analytic and numerical results for a 
plasmon being accelerated th rough the interaction with a wind 
(see e.g. jNittman et al. 1982[ ). 
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Fig. 4. Density and velocity at t = 750 years. See Fig. || for more details 
Table 1. Line strenghts (in Lq) 



10 



17 



Time (yrs) 


Ha 


[N II]A6583 


[N II]A5755 


[S IIJA6731 


[S IDA6717 


250 


1.186 x 10" 3 


5.310 x 10" 4 


8.691 x 10" 6 


1.224 x 10" 4 


8.799 x 10" 5 


500 


1.025 x 10 -3 


1.581 x 10 -4 


3.069 x 10~ 6 


8.937 x 10" 5 


6.541 x 10~ 5 


750 


8.896 x 10~ 4 


7.773 x 10~ 5 


1.502 x 10" 6 


4.252 x 10~ 5 


3.145 x 10~ 4 



Table 2. Peak line ratios to Ha 



Time (yrs) 


[N II]A6583 


[N II]A5755 


[S II]A6731 


[S II]A6717 


250 


2.000 


0.0442 


0.636 


0.454 


500 


1.498 


0.0308 


0.844 


0.559 


750 


2.047 


0.0436 


1.112 


0.770 



6. Comparison with observations 

The number of ions included in our numerical simulation 
makes it relatively straightforward to produce synthesized ob- 
servations of various kinds which can be compared to the ob- 
served FLIERs. In Table 1 we show some line strengths derived 
from the snapshots at t = 250, 500, and 750 years. The top 
part lists the total luminosities (in solar luminosities) emitted 
by the clump material. We limit ourselves to the clump mate- 
rial to exclude the large amount of Ha emission produced by 
the environment. The lower part of the table shows the maxi- 
mum ratio of two forbidden lines of [N II] and [S II] to Ha, 
derived using the whole image, in this case including the envi- 



ronmental Ha emission. However, the resulting ratios are not 
very different if one excludes the contribution from the environ- 
ment. The model reproduces the main observational character- 
istic of FLIERs and ansae, namely that they are very bright in 
[N II]A6583. Note that we use 'normal' nitrog en abundances 
in t he model and therefore supp ort the view of Dopita (1997) 
and [Alexander & Balick (1997)| that the high [N II]/Ha ratios 
do not necessarily imply high nitrogen abundances, but can be 
produced in these type of ionization fronts. 

The other line strengths are also not far off from the ones 



found by Balick et al. (1994). Using the ratio of the red [S II] 



lines one can derive an average electron density, which turns 
out to be around 1500 cm~ 3 at all epochs. Similarly the ratio 
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Fig. 6. Left: Contour plots of the pressure, temperature and H ionization fraction at t = 500 years. The dashed line shows the fit 
of Eq. [28] with h = 2.7 x 10 15 cm. The contours in the pressure and temperature plots are logarithmic with a factor \/2 between 
each contour and minima of 1 x 10~ 10 erg cm" 3 and 100 K. The contours in the ionization plot are linear from 0.02 to 0.88 in 
steps of 0.02. Right: Cuts along r = for the three contour plots. The dashed line in the pressure plot shows the pressure gradient 
corresponding to the scale height found from fitting the shape of the plasmon. 



of the two [N II] lines gives an electron temperature, which is 
around 10 4 K. The density value matches the observed one very 
well, but the temperature appears to be slightly higher than the 



observed one (cf. Balick et al. 1994) 



Figure [7] shows images of the plasmons in the Ha and 
[N II]A6583 lines projected at an angle of 0° in the sky. The 
typical morphology is that of a little arc pointing away from 
the star, both in Ha and [N II]. However, in Ha the contrast 
between the plasmon and the environment is much lower than 
in the [N II] line, making the object less noticeable, which is 
in accord with the observations of FLIERs and the other low- 
ionization knots. [O III]A5007 (not shown) is similar to Ha, 
again as observed. The morphology matches less well. The 
model images show arcs pointing away from the star, whereas 
in the HST observations the FLIERs appear to have mostly 
arcs pointing towards the star, although there are also some in- 
stances for which the reverse appears to be the case. Also, the 
cones pointing back towards the star as seen in NGC 6826 and 
NGC 7009, have a reverse orientation compared to the cones 
seen in our Ha images. However, the HST images also show 



clearly that the real FLIERs are not homogeneous knots such 
as in our model, making the comparison difficult. 

Finally we show in Fig. || the profiles of the [N II]A6583 
line for the three epochs that we are considering. The instru- 
mental broadening is assumed to be 7 km s _1 and the objects 
are moving away from us at an angle of 45°. At t = 250 years 
the instrumental broadening is enough to hide any motion 
of the plasmon, but at the two later times one sees how the 
line becomes shifted due to the motion of the plasmon. At 
t = 500 years this shift is about 7 km s _1 , at t = 750 years, 
17 km s _1 . At t = 750 years the line profile clearly deviates 
from thermal, showing a low velocity tail. The observed line 
profiles in [N II] are typically 30 km s _1 wide, perhaps due to 
the fact that we see an ensemble of evaporating clouds. 

Assuming that the clump is initially moving with the same 
velocity as the environment, this environmental velocity should 
be added to obtain the space velocities of the clump. The value 
of the environmental velocity depends on the kinematics of the 
nebula, but for elliptical PNe can be estimated to lie between 10 
and 50 km s _1 , which gives absolute FLIER velocities in the 
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Fig. 7. Ha and [N II] images for times 250, 500 and 750 years. 
The contours span a dynamic range of 1000 in steps of a factor 
2 



range 30-70 km s _1 , close to the observed values. As we saw in 
Sect. 5, the average velocity of the plasmon is about 26 km s _1 . 
This average velocity is achieved when 2.8 x 10 -6 M© (12.5% 
of the original mass) is left. Higher velocities are achieved later 
on, but by this time only a tiny fraction of the original mass 
is left. The plasmon reaches a velocity of 100 km s _1 when 
only 1% of the original mass remains. One would need to start 
with very massive clumps (~ 10 -4 Mq) in order to explain 
FLIERs whose velocity differs by more than 30 km s _1 from 
their immediate environment. For those an additional acceler- 
ation mechanism is needed, for example a high initial clump 
velocity or a fast wind flowing past the clump. We intend to 
explore models of these types of flow in future papers. 

7. Conclusions 

We studied the evolution of a dense clump in the radiation field 
of a hot central star of a PN using both a simple analytic ap- 
proach and a detailed numerical model. We compared the re- 
sults to each other and to the observations of FLIERs. The main 
conclusions from this are: 

1 . The evolution of an accelerating, photo-evaporating spher- 
ical clump can be described by a relatively simple analytic 
model which gives both the velocity, the mass and the shape 
of the clump with time. 




-30 30 

v r (krn/s) 

Fig. 8. Line profiles for [N II] A6583 at times 250, 500 and 750 
years. The inclination of the flow direction with respect to the 
observer is taken as 45° and the instrumental broadening is 

7 km s _1 



2. When comparing this solution with a detailed numerical 
model we find an excellent agreement, implying that the 
simple analytic model does indeed give a good description 
of the evolution of such a photo-evaporating clump. 

3. The synthesized observations compiled from the numeri- 
cal model can reproduce many features of FLIERs, such as 
typical line strengths (especially the strong [N II] lines), 
the average velocity difference with the environment, and 
perhaps the morphology. The reproduction of strong [N II] 
lines is noteworthy since we use standard abundances. We 
therefore conclude that the photo-evaporating clump model 
is moderately successful in explaining the FLIER phe- 
nomenon. However, some of the morphological features 
seen in the HST observations definitely need further study- 
ing and modelling. 

The model we studied here does not address the origin of 
the clump. Given that the FLIERs lie along the major axis 
of the PN, it is logical to assume that their origin is related 
to the cause of the asphericity of the nebula. There are then 
two options, either the clump was formed very near the star, 
perhaps through some process involving a disk or a jet (see 
|Soker 1996 and references therein), or during the early phases 
of the interaction between the stellar fast wind and the sur- 



rounding material (Frank et al. 1996). The first explanation re 



ceived some support from the reported high nitrogen abun- 



dances, but as was shown here, as well as in Dopita (1997) and 



Alexander & Balick (1997), the nitrogen abundance does not 



need to be exceptional if one assumes the FLIERs to contain 
an ionization front. Recent NICMOS results on the proto-PN 
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CRL 2688 seem to indicate that dense knots indeed form dur- 
ing the proto-PN phased]. 

The analytic solution we derived applies to any photo- 
evaporating clump, and could therefore be applied to other 
cases such as the cometary globules found in the Helix neb- 



ula (O'Dell & Handron 1996; Meaburn et al. 1997 1, and dense 



knots near massive stars in star formation regions. 
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